#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018/1/31 17:00
# @Author  : Aries
# @Site    : 
# @File    : interplateImage.py
# @Software: PyCharm Community Edition

from image_tools import *
from config import *
import numpy as np
from scipy.spatial import cKDTree
from osgeo import gdal



def cKDTreeLookup(imagefile, demfile, lutfile, outdir):
    """
    :rtype is no return
    """
    # read radiance image and dem tiff
    padfTrans, b2, b3, b4, b5, b6, b7= readRadianTIFF(imagefile)
    #print('minx, maxy:', padfTrans)
    print('minx, maxy:', padfTrans)
    print('image size:', b2.shape)
    #print(b2[1500, 2000])
    dem = readDEM(demfile)
    dem = dem[:-1, :-1] #same size
    print('dem size:', dem.shape)

    # stack coords line by line
    ## so we need to interplote line by line
    ### this is also efficent to save array as image!
    coords2 = np.stack((b2, dem), axis=-1)
    print('stacked size:', coords2.shape)
    coords3 = np.stack((b3, dem), axis=-1)
    coords4 = np.stack((b4, dem), axis=-1)
    coords5 = np.stack((b5, dem), axis=-1)
    coords6 = np.stack((b6, dem), axis=-1)
    coords7 = np.stack((b7, dem), axis=-1)


    #load look up table
    lut_table = np.genfromtxt(lutfile, delimiter=',', skip_header=1, dtype='f')
    print('size of look up table:' ,lut_table.shape)

    #6S computed surface reflectance for each bands
    s2 = lut_table[:, 2]
    s3 = lut_table[:, 3]
    s4 = lut_table[:, 4]
    s5 = lut_table[:, 5]
    s6 = lut_table[:, 6]
    s7 = lut_table[:, 7]

    #load LUT dimensions
    LUT = lut_table[:, 0:2]
    print('shape of kdtree dimensions:', LUT.shape)

    # build Tree based on LUT fields
    tree = cKDTree(LUT)


    #query tree and interpolate for band2
    d_b2, i_b2 = tree.query(coords2, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b2].T
    #interpolation
    v_b2 = ((v_k2 + v_k1) / 2).T
    print("result image array shape:", v_b2.shape)
    #print(v_b2[1500, 2000])
    array2raster('{}/LC08_L1TP_150033_b2_surf.tif'.format(outdir), padfTrans, v_b2)  # convert array to raster

    #query tree and interpolate for band2
    d_b3, i_b3 = tree.query(coords3, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b3].T
    #interpolation
    v_b3 = (v_k2 + v_k1) / 2
    array2raster('{}/LC08_L1TP_150033_b3_surf.tif'.format(outdir), padfTrans, v_b3.T)  # convert array to raster

    #query tree and interpolate for band2
    d_b4, i_b4 = tree.query(coords4, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b4].T
    #interpolation
    v_b4 = (v_k2 + v_k1) / 2
    array2raster('{}/LC08_L1TP_150033_b4_surf.tif'.format(outdir), padfTrans, v_b4.T)  # convert array to raster

    #query tree and interpolate for band2
    d_b5, i_b5 = tree.query(coords5, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b5].T
    #interpolation
    v_b5 = (v_k2 + v_k1) / 2
    array2raster('{}/LC08_L1TP_150033_b5_surf.tif'.format(outdir), padfTrans, v_b5.T)  # convert array to raster

    #query tree and interpolate for band2
    d_b6, i_b6 = tree.query(coords6, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b6].T
    #interpolation
    v_b6 = (v_k2 + v_k1) / 2
    array2raster('{}/LC08_L1TP_150033_b6_surf.tif'.format(outdir), padfTrans, v_b6.T)  # convert array to raster

    #query tree and interpolate for band2
    d_b7, i_b7 = tree.query(coords7, k=2)
    #find 2 nearest neighborhoods
    v_k1, v_k2 = s2[i_b7].T
    #interpolation
    v_b7 = (v_k2 + v_k1) / 2
    array2raster('{}/LC08_L1TP_150033_b7_surf.tif'.format(outdir), padfTrans, v_b7.T)  # convert array to raster





def interpolate():
    """
    interpolatation
    :param inpath:
    :return:
    """
    # readfile and call interpolate function
    inpath = "/home/deyu1/share/gazreasarch/optical_related_experiments/LandsatData/6Satmos/"
    imagefile = inpath + 'LC08_20160227_TOA_6s_clipgeom.tif'
    demfile = inpath + 'Fill_clip_Gaiz_toUTM_30m_snaped.tif'
    lutfile = inpath + 'LC08_L1TP_150033_LUT_BRDF.csv'
    outpath = "{}/LC08_20160227_surf".format(inpath)
    cKDTreeLookup(imagefile, demfile, lutfile, outpath)



if __name__ == '__main__':
    interpolate()